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ABSTRACT 


An algorithm is developed which enables the computation 
of horizontal range and/or depth from a submerged sound 
Source, usimigq ray acomst@cs and the Lloyd's mirror 
interference effect. The solution is based on Snell's law 
and involves integrating multipath sound rays to find the 
difference in length between the direct and surface 
reflected sound paths from the source to the receiver. This 
difference in path length is directly related to the 
observed Lloyd's mirror interference pattern. 

No assumptionsas to the mathematical characteristics of 
the sound velocity profile (SVP) are made nor are far field 
approximations used. The solution is as accurate as the SVP 
data input to the problem. 

A general computer flow chart and basic language program 
are provided to allow local commands the capability of pre- 
mission planning based on specific operating area 


environmental information. 
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I. INTRODUCTION 


Sound propagation between a source and a receiver in the 
ocean may occur by many different ray paths, each dependent 
on existing physical and acoustical conditions. Crucial to 
the determination of which sound paths are of importance in 
any given problem are factors such as the sound velocity 
profile, the depth of the water, the bottom contour and 
composition, the sea surface conditions and the intensity 
and frequency of the sound. 

This thesis deals with two specific paths of sound 
propagation from a source to a receiver. The first, called 
the direct path, does not touch the surface or bottom of the 
sea as it travels from the source to the receiver. The 
second path reflects off the sea surface before reaching the 
receiver and is called the reflected path. 

The following explanation of the reception of sound from 
a Single source via these two distinct paths is taken 
directly from Urick: [Ref. 1: pp. 120-123] 

"Tf the sea surface were perfectly smooth, it would 
form an almost perfect reflector of sound. The intensity 
of sound reflected from the smooth sea surface would be 
very nearly equal to that incident upon it." 

He continues: 

‘ee ceiterion forthe roughness or smoothness of a 

surface is given by the Rayleigh parameter, defined as 


R= kH sin@ ;where k is the wave number, (27)/\A (A is the 
wave length of the sound incident on the surface), H is 


Dudl 
Mont 


the rms (sea) wave height (crest to trough), and 8 is the 
grazing angle”. 


Sea Surface 





R = kH sin 8 


Brgure I-11, Rayleigh Parameter (R) for Describing Sea 
Surface Roughness. 


"When R<<1, the surface is primarily a reflector and 
produces a coherent reflection at the specular angle equal 
to the angle of incidence. When R>>1, the surface acts 
as a scatterer, sending incoherent energy in all 
Gacect ions. When certain theoretical assumptions are 
made, (See: P. Beckmann, and A. Spizzichino: Scattering of 
Electromagnetic Waves From Rough Surfaces, p. 93, The 
Macmillan Company, N.Y., 1963.) the (pressure) amplitude 
reflection coefficient w of an irregular surface (defined 
as the ratio of the reflected or coherent pressure 
amplitude of the return to the incoherent pressure 
amplitude) can be shown to be given by the simple 
expression w= exp (-R)." 


Urick goes on to explain the image interference. 


"When the sea surface is not too rough, it creates an 
interference pattern in the underwater sound field. fThis 
pattern is caused by constructive and destructive 
interference between the direct and surface reflected 
sound and is called the Lloyd mirror, or image 
interference, effect.” 


In Urick's sound propagation in the sea [Ref. 2: p. 9-5], he 
continues the discussion of surface reflection interference: 


"This sound field may be divided into three parts, 1) a 
NEAR FIELD close to the source in which the image source 
is too far away, and the reflected sound is too weak to 
produce appreciable interference, 2) an INTERFERENCE FIELD 
in which there are strong loops and nulls in the signal 
received by a receiver moving outward in range, and 3) a 
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FAR FIELD in which there is an increasingly out-of-phase 
condition between source and image and the intensity falls 
off as the inverse fourth power of the range (transmission 
loss increases as 40 log r)." 
The interference effects felt at the receiver is 
expressable as the ratio of the combined direct and 
reflected sound intensities to the direct sound intensity 


alone. 


Stated mathematically: 


a ee . 
tal d i 
Intensity Ratio = tota Deedneneciumene! Tecticd) 
“Gli eeke Tdirect 
or more precisely 
2 Rp\ Rp 
Pmeensity Ratio = T+ up —_ =o) ey COS (Wy) 
Rp Rp 


Where Rp is the direct ray path length, Rp is the reflected 
path length, uw is the surface pressure reflection 
coefficient, w (where w = 2rf) is the angular frequency of 
the sound, and +r is the time delay caused by the difference 
in the two path lengths. A formal derivation of the 
intensity ratio is given in Appendix A. A more 
comprehensive treatment of the effects of surface roughness 


on acoustic waves is given by Clay and Medwin [Ref. 3]. 


V1 
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II. APPLICATION of LLOYD MIRROR INTERFERENCE PHENOMENON to 


SOURCE FROM the RECEIVER 


Tactical tracking information from the Lloyd's mirror 
phenomenon has been used operationally by the Navy for many 
years. Inan excellent report on the analysis of Lloyd's 
Mirror phenomenon under constant velocity gradient 
conditions, Shudde [Ref. 4: p. 3], discusses the difficulty 
of solving the surface reflection interference problem. 
Accurate error estimation for any acoustical condition other 
than the isovelocity case is not practical unless a large 
computer is used. Since Shudde's 1977 report, the 
advancements in memory size, and speed of operation of small 
portable desk top computers now allow an accurate solution 
of source position for all sound velocity conditions. The 
following mathematical solution is based on the lengths of 
the sound paths determined from Snell's Law. As explained 
byewewek. (Ref. 1: p. 116], 

"One of the most important practical results of ray 
theory is Snell's Law, which describes the refraction of 
sound rays in a medium of variable velocity." 

The solution is not dependent on any isovelocity or constant 
gradient conditions other than the reasonable assumption 
that no horizontal sound speed gradients exist throughout 
the operating area. The mathematical solution presented in 


Section III makes no assumption as to the mathematical 
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Cmaemacteristics of the sound velocity profile nor does it 
use far field approximations. Past solutions depended on 
specific sound velocity profiles, such as an isovelocity or 
enact laine’ protale wwhich resulted in) significant 
positional errors for some tactical applications when the 
actual sound velocity profile did not match the predicted 
model. 

PieaPpivangesthne solution to the Lloyd's mirror 
interference effect one must understand the relationship of 
the propagating sound to the interference pattern. Coppens, 
Sanders and Dahl [Ref. 5 : p. 55] explained, 


"...the noise radiated by underwater vehicles consists 
of two contributions: broadband noise (Similar to the 
hissing found between FM radio stations) and tonals 
(single frequency tones, like whistles). Thus, the 
radiated noise consists of all frequencies between certain 
broad limits with certain specific frequencies occurring 
with high intensity. Each of these frequencies, whether 
of broadband noise or tonal, has its own interference 
pattern. 

Because each frequency has its own period and 
wavelength, each will have a different phase delay for 
propagation over the direct and surface reflected paths. 
As a result, the ranges between source and receiver for 
which the surface reflected signal tends to cancel the 
direct signal will be different for each different 
frequency. Thus, the ranges at which different 
frequencies will have nulls at the receiver will be 
different. This means that as the range between source 
and receiver changes, the frequencies which have nulls at 
the receiver will also change." 


Figure II-1 shows a schematic time vs. frequency plot of the 
interference pattern as the distance between the source and 
receiver is decreased and then increased. inecdistanct 


tonals will most readily show the effects of Doppler shift 
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(apparent frequency change) as minimum separation range, is 
passed (see t = 0 in Figure II-1). The light (unshaded) 
hyperbolic shaped regions centered about t = 0 are known as 
the Lloyd's mirror interference pattern for the broad band 
noise. The frequency values Fy, F5, F3, Fy, represent the 

Minimum frequency of each hyperbola. The actual frequency 
from the sound source that is being nullified by the surface 
reflection interface varies with time as the hyperbolas are 
traced out. The measured separation in hertz between any 
two successive nulls at any time t, is termed the difference 
frequency, Af, and is the frequency value required in the 
problem solution discussed in the following Section. (See 


Section III.D.2) 
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Figure II-1. Frequency vs. Time Display of Lloyd's Mirror 
interference Pattern 
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III. RANGE vs. DEPTH SOLUTION USING the LLOYD'S MIRROR 


INTERFERENCE PHENOMENON 


A. THEORETICAL BACKGROUND 

1. Snell's Law 

Snell's Daw, (Ref. 1: p. 116, Ref. 63° p. 120], is a 

result of Fermat's minimum time principle which states that 
each ray will travel the path requiring minimum time of 
travel from source to receiver. 

Mathematically, Snell's Law is expressed as: 

sin 9 


Snell Constant = ——— : G) 
ez , 


Where § is the angle that the ray makes with the vertical 
and c(z) is the sound velocity for a given depth z. (See 


Figure III-1) 





aljglas ie ales ill 
Surface + Range, 
Meters 
Raypath direction at Depth 4@,, 
in oe a as al 
Depth (2), 
Meters 





Figure III-1. Snell's Law Applied to Initial Sound Source 
Ray Angle. 
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The sound ray maintains the same Snell Constant value 
throughout its transmission path. Sound velocity changes 


result in an increasing or decreasing angle 9. 


2. Ray Path Integration 


Expanding the diagram of the sound ray, the 
differential raypath length, ds, can be described as a 
function of the differential change in depth, dz, and the 
differential change in horizontal range, dr. (Figure III-2) 

Sound 
Ray 


ar / 
dz ds ds = dz? ~ dr? 
Source 


Figure III-d2. Pythagorean Theorem Applied to Ray Path 
Length. 

To calculate the actual path length of the sound 
ray, smaller and smaller segments of path length, As, 
horizontal range, Ar, and depth, Az are used until the 
differential segments ds, dr, and dz are approached. This 
tends to smooth the curve of the ray into a sequence of 
straight line segments. The path length traveled by a 
particular ray can then be determined by integrating between 
the source depth and receiver depth. The possibility of the 


direct ray crossing the same depth more than once is not 


ile 





Somstaercd a problem in this solution due to the short 
ranges involved in Lloyd's mirror observations and the long 
ranges required for a sound ray to curve back through the 
same depth before reflecting off the surface. To handle the 
fact that the surface reflected ray does cross the same 
depth twice in the ranges pertinent to this discussion and 
to avoid errors in integration, the concept of a geometric 
equivalent image [Ref. 1: p. 123], of the source is used. 
The image reflected path and the actual reflected path are 
the same length. (Figure III-3) 


' 


Bj; Image Reflected Path 


“Zp ~ 
AN 
‘XN 
~ 
ae & 
ez) Surface o a +Range © 
SVP Meters 
Source 
+Z), Actual 
B Reflected 
Path 
+Z, A 
+Depth 
Meters 


Figure III-3. Actual and Image Sound Paths From Source to 
Receiver for a Given Sound Velocity Profile 

The curved sound path shown in Figures III-3 and 

ITII-4 is due to the change (in general) in the sound 


velocity as the depth changes. It is possible to show that 


Enewpatn £rom A to B is the same’length as the path from A 
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EOmo Also, bY inspection of the sketch, the horizontal 
range from A to B or A to B. 1s the same. Therefore, for 
the reflected path calculations of path length and 
horizontal range, the image reflection depth of aA dM [sXe 


used. 


B. HORIZONTAL RANGE SOLUTION 
(Pete ct Pati 
The horizontal range traversed by the direct path 
sound ray can be calculated from Snell's Law and the 
Pythagorean Theorem for right triangles applied to Figure 


III-2. 
From (1) , the Snell Constant, K, is: 


sin @,; 


K, = —____— 


c(z;) 


This is the Snell constant for the particular direct path 
ray that ultimately travels from the source to the receiver, 


but initially starts at depth z; and moves in direction 4,. 


Using the Pythagorean Theorem on Figure III-2 


dr 


sin @ = (2) 


\/ (dr)* + (dz)? 


Solving for sin @, where @ = variable for the 
particular direct path ray 
tagged with Snell constant 
K 

17 
sin @ = K, c(z) (3) 





he, 








setting (2) and ©) equal to each other yields: 


ar 


Ky c(z) = (4) 


\/ (dr)* + (dz)? 


The right side of (4) can be simplified by dividing 


numerator and denominator by dr. 


Ky Cuz) = 
1 + (=) 
ar 


Solving for dz/dr follows in steps: 


dz 2 1 
ips -_ pees cree 
dr K, ¢(z) 


dz\ ¢ 1 

1 ge (ee > —_—_—, 
Gi [Ky ez) 
(as) 1 ~ (K, ¢(z)]? 
ar : [Ky c(z)1? 
daz 
= é 
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dz 
dr 
Now, solving 
range yields: 
aig 


= Ky ez) 


fourdr, 


Ky c(z) 


/i - (Ky c(z)]¢ 


+ i he eC 
We [hse eli ic 


© 


the incremental horizontal 


az 


©) 


Integrating (6) with respect to depth, with the 


limits of integration being the depth of the source and the 


depth of the receiver, 


Rea adirect path) = 


The “Range" 





yields: 
source 
Depth 
Ky c(z) 

+ 

Se ra4)) 
Actual 
Receiver 
Depth 
(along path of ray 
"tagged" K, ) 
is always positive 


dz 


and thus the 


appropriate + or - sign is used to make this integral always 


positive regardless of the integration limits. For ease of 


discussion in the rest of this solution, the receiver will 


be assumed to be at a shallower depth. 


Za 





2. Reflected Path 


Recalling the image depth discussion and noting a 
See cent eaagre, ~O., for the reflected ray is required 
(Figure III-4), which results in a different Snell constant, 


K5, where 


Sin (6; ) 
c(Z; ) 
It follows directly that (7) can be rewritten for the 


reflected horizontal range as; 


Source 
Depth 


Ky c(z) 
2 
R (reflected path) = — —————————_- dz (9) 


Vibes x ale) 4 


Image 
Receiver 
Depth 


ZZ 





R( reflected Sey ee 
| 


as Image Receiver | 
oe | 


at + Range 
O - Meters 
24 , sin (0) 
K, = ———— 
+Zy, Receiver | c(z.) 
} | Sin (8) 
| K> = 
4 Sz} 
Depth sill Ze : 
Meters R(direct path) ——————4 


+ 
Figure III-4. Common Direct and Reflected Sound Paths From 
Source to Receiver. 
C. SOUND RAY PATH LENGTH SOLUTION 
lTeberect Path 
The direct path ray length can be calculated by 
Beeeehins to Figure III-2, applying the Pythagorean Theorem 


and noting: 


ds = \/(dr)2 + (dz)? 


realizing the path length S from receiver to source is given 
by: 


Source Source 
Depth Depth 


S = ds 
| | \| (dr)? + (dz)? a) 


Receiver Receiver 
Depth Depth 
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dz 





Which can be simplified by multiplying by 


Ou 
N 
NO 


as follows: 


Source 
Depth 
S (Direct Path) = 1+ /ar\ 2 
oe. Gz 
Vaz 
Receiver 


Depth 


dr 


Substituting the equivalent of (=) from the reciprocal of 
dz 


©, 


Source 
Depth 


S (Direct Path) = 
dz 





Receiver 
Depth 


Source 
Depth 


S (Direct Path) = 
diz, 





Receiver 
Depth 
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source 
Depth 


S (Direct Path) = —_-- - dz 


1 - [Ky c(z)]* 


Actual 
Receiver 
Depth 


2. Reflected Path 


Using the same arguments about the reflected path as 


presented for horizontal range and substituting K>5 for Ky 


and image receiver depth into (17) yields: 


Source 
Depth 


S (Reflected Path) = ———— dz 
Wi 28fk, e@z)12 


Image 
Receiver 
Depth 

D. EXPLANATION OF SOLUTION 


1. Some Observations 


Equations (7) and (12) are path integrals for the 


direct ray, using the same Snell constant and integration 


limits. Further, equations (9) and 03) are path integrals 


of the reflected ray and thus use a different Snell constant 


from the direct path integrals, but use the same integration 


Dns 








limits respectively for the horizontal range and ray path 
lengths. 
Ze) SOluUtion 

The key to the eventual numerical solution is the 
equality between the direct path horizontal range and the 
reflected path horizontal range. Although the actual path 
length of the reflected ray is longer than the direct path 
ray between common source and receiver, the horizontal 
distance traveled by both is the same. 

By calculating all of the initial ray angles that 
are possible for each depth of the source of interest, two 
horizontal range arrays, one for the direct path and one for 
the reflected path, are formed. Each of these arrays is 
linked respectively to a direct and reflected ray path 
length array by the corresponding Snell Constant used to 
determine both the horizontal range and sound ray length for 
each path, (See Figure III-4). For example, suppose the 
initial direct path ray angle of 15 degrees from the 
vertical covered a horizontal range of 20 meters in 
traveling from the source to the receiver. Similarly, the 
initial reflected path ray angle of 8 degrees from the 
vertical covered the same horizontal range of 20 meters from 
the source to the receiver. Then the difference between the 
reflected ray path length for 8 degrees and the direct ray 


path length for 15 degrees could be computed (Figure III-5). 
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20 meters 





+ 
Range 
Meters 


Zh 
8 degrees 
15 degrees 
Zs Source 
+ Depth 
Meters 


Figure III-5. Example of Reflected Path Ray Interfering 
With Direct Path Ray at the Receiver. 


ir @Sethaissdifftference in ray path length, 


AS = Sp = Spe 


time delayed phase shift sensed at the receiver and known as 


for a wide band of frequencies that causes the 


the Lloyd's mirror interference phenomenon. Since both the 
direct and reflected rays are traveling at the same sound 
Voilocmevmateche receiver, the ray path length difference can 
be directly converted to the frequency difference, Af as 
observed at the receiver and seen in the Lloyd's mirror 
interference pattern. 


The intensity ratio at the receiver is given by 


Rn \* R 
2 D D 
Intensity ratio = ] +n — - 2ui — cos (wT) 
Rr Ss 
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The ratio reaches a minimum value corresponding to 
maximum interference when cos(wr) = 1. This occurs’ when 


wr= 2" or any integer multiple of 2nr. 


Note that = c/f where \X = wavelength (m) 
c = phase speed (m/s) 
f = frequency in (hz) 
k = wave number (m~!) 


syovel US 7A ON 

Now let 4% equal the phase difference seen at the 
receiver due to the two path length differences and defined 
as, 

A®G=zk (Sp- Sp), where 

Sp = Path length of the reflected ray and 

Sp = Path length of the direct ray. 

Let AS = (Sp- Sp) 

So that the first null occurs when 

Bape & (45) = 27 = wr 

So k (AS) =wr= 2wr=27(A4S)/X or 

AS = 2rd\/2*” = 

therefore AS = c/f, 

so the first null occurs at frequency 


f, = c/ AS 


The second null occurs at 
wr=- 47 so that 
AS = 4rrA/2r= 2X = 2c/f5 


therefore second null occurs at frequency 
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and 


where m >> n for best accuracy so that the measured 
difference in frequency between any two successive nulls can 


be related to the difference in path length required to 


£5 = 2c/AS = 2(f,) 
which shows that £5 _ f, = f, 


and so on. So that for the nth null 


fe = nce/AS = nf, 
fi - fp-1 = £1, or more generally 
dees = fs C 
m n 
= f, = lo. 
m-n 42S 


produce that null by the equation 


C(receiver ) 


AS 


f= 


where Ci receiver) is the actual sound velocity at 


the receiver at the time of reception. 


28, 





IV. COMPARISON of the ANALYTICAL and COMPUTER SOLUTIONS for 


Se a Rt em ee 


This analysis demonstrates the accuracy of the computer 
solution by using an isovelocity sound propagation model. 
The constant sound velocity condition results in straight 
line transmission paths as the sound travels from the source 
to the receiver. By solving the geometric model analytically 
and comparing the results to the computer output one can see 
the accuracy of the Romberg integration technique coupled 
with the four integrals derived in section III . The 
isovelocity case is the simplest propagation model yet it 
exercises all phases of the computer program's main 
calculation algorithm. The only section of the overall 
program not tested by the isovelocity case are the sound 
velocity calculations performed as the user inputs a varying 
temperature profile. Table IV-1 compares the difference in 
path lengths between the direct and reflected sound paths as 
calculated analytically and by the computer. Of the 20 data 
points compared, the computer missed the exact ray path 
length for a horizontal range of 750 meters by only 19.16 
centimeters in the worst case and averaged 5.02 centimeters 
error over all 20 cases. Table IV-2 shows the resulting 


horizontal range error due to the error of the computer's 
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> Gane Eft) 


calculation presented in Table IV-1. The average range error 
for the 20 points is only 1.52 meters. 

The computer program steps through sea depths of 0 to 
200 meters in 10 meter increments where each depth 
represents a possible target depth. To calculate a complete 
set of Snell constants, initial sound ray angles of 3 thru 
177 degrees in increments of 3 degrees ( measured from the 
vertical ) are used. Ultimately the integration section 
derives a path length difference between the surface 
reflected ray and the direct path ray for each of these 
selected starting ray directions.The 59 initial rays times 
the 20 depth increments yields a total of 1,180 data points. 
These point solutions are used to interpolate values of path 
length difference for points every 10 meters of horizontal 
range out to 1000 meters which results in a two dimensional 
matrix containing 2000 ( 20 depth increments times 100 range 
increments ) horizontal range vs depth elements. 

To compare the computer solution to actual geometric 
results, five ranges for each of four depths were chosen at 
random. In all comparisons the receiver is set at a depth of 


27 meters (approximately 90 feet). 
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SOURCE 
DEPTHS 
METERS 


10 


50 


90 


e.0 


170 


TABLE IV-1 


ISOVELOCITY CASE 





COMPARISON BETWEEN 


GEOMETRICAND COMPUTER SOLUTIONS 


lovee) 0 7 


15282728 


AS01755 


Aro 196 


BL 5590 / 


S20 6 


5) Gre Gee ye) 


53.3808 


53.6205 


53.9714 


PATH LENGTH DIFFERENCE 


ANALYTICAL 


COMPUTER 


HORIZONTAL RANGES, METERS 


100 


5.903 


2 pnd Fee OB, 


2395992 


23262020 


S572 10 


BSechas 


42.5840 


42.6136 


46.4312 


46.4396 


380 


1. 4170 


1.4418 


7.0274 


ieV3SC9 


V2 47 


12.4688 


17.4439 


17.5644 


DEX We 72 


22.0264 


71590 


Ona 915 


(Ses &: 


Bmp 9 7 


3.5714 


6.4289 


G26 7 


A ss, 


93.0255 


Alte a0 2 


Va oo 





TABLE IV-2 


COMPARISON of PROGRAM RANGES vs. ACTUAL RANGES 


ACTUAL 
SOURCE 
DEPTHS BORIZONTAL RANGES, METERS 
METERS 
20 100 380 750 
10 
19.954 99.935 B73. 42 S06. 3 
50 
20.134 Soo lente fis OW os.3 7 7] S53 nor 
90 
mon. 745 100.038 SiGe oS 11a Oe 
VS9) 
19.42 99.815 S77 e075 766.41 
170 
UA PAS 99.93 Sel co 760.475 
COMPUTER 


HORIZONTAL RANGES, METERS 
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V. COMPUTER PROGRAM FLOW CHART 


A. FLOW CHART and PROGRAM BACKGROUND 

The following block diagram flow chart demonstrates the 
specific steps required to obtain an accurate horizontal 
range vs. depth solution using the Lloyd mirror interference 
phenomenon. Specific values or data used in the included 
program will be given at the bottom of each appropriate 
block. 

The actual numerical integration section of the program 
is taken directly from the Romberg integration program 
provided by Miller [Ref. 7: pp. 215-218]. 

The sound velocity equation used in the program was 
developed by Coppens [Refs. 6 and 8] and selected after 
several equations of varying complexity were considered. 
Specifically, equations derived by Wilson and Leroy [Ref. 1: 
p. 105], plus Mackenzie [Ref. 9] were compared with Coppens. 
All results were within 1 to 2 centimeters per second of 
each other. Coppen's equation was chosen because of the 
ability to correct for not only temperature, pressure and 


salinity, but also, latitude variations. 
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B. FLOW CHART 
1. Initialization and Main Menu 


INITIALIZATION 
DEFINE ALL VARIABLES 
2. DIMENSION ALL ARRAYS 
3. DEFINE TOLERANCE 

DESIRED FOR NUMERICAL 
INTEGRATION 


yo | 
e 


RANDOM ACCESS FILE 
CREATE A RANDOM 
ACCESS FILE. TC 
STORE THE RANGE 
Won DEPTH BMATR Ix 
SOLUTION 





RETURN FROM RETURN TO 
ERROR STATEMENT MENU FROM 
FOLLOWING PROGRAM 
INTEGRATION ERROR SUBROUTINES 
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DISPLAY MENU OF 
PROGRAM SUBROUTINES 
ALLOWING USER TO 
CHOOSE DESIRED 
OPERATION AS 
FOLLOWS: 
REINITIALIZE RANDOM 
ACCESS ANSWER FILE 

AFTER INITIAL 

RUN AND PRIOR TO 
ALL SUBSEQUENT RUNS 
INPUT WATER TEMPERATURE 
INFO AS A FUNCTION OF 





















3. INPUT SEA SURFACE TEMP- 


ERATURE AND AVERAGE ====15) 


GRADIENT 


4. INPUT SOUND VELOCITY 


PROFILE AS A FUNCTION +4) 


OF DEPTH 
[ 5. COMPUTE HORIZONTAL RANGE ees 
VS. DEPTH SOLUTION AND 3) 
STORE IN ANSWER FILE 
6 eb eSceibAY SANSWER FILS ON 
nae Menno Creve OR lo. 
ON A PRINTER 


EXIT PROGRAM AND RETURN 
DO OPERATING SYSTEM 
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2. Reinitialize Random Access Answer File 


——— ES SS 





WARN USER THAT 


THIS SECTION WILL | RETURN 
DESTROY CURRENT ' NO TO MENU 
ANSWER FILE 
QUERY USER 


REINITIALIZATION 
DESIRED, YES OR NO 









Bele ont oe Fie 
GCLOSOMANSWER FILE 


KIEL @B=ECETE) ANSWER i ecl 
OPEN NEW ANSWER FILE 





RETURN 
TO MENU 





oy 





3. Input Water Temperature. Profile 


© 


INPUT SALINITY IN PARTS 
PER THOUSAND FOR OPERATING 
AREA, ENTER A VALUE OF 35 

IF SALINITY IS UNKNOWN 














INPUT LATITUDE IN WHOLE 
DEGREES FOR OPERATING AREA 


Pee Wao Ce niPe nA PUR Geeks 
DEGREES CELSIUS AT UNIFORM 
| Nenwretwro OV Bel TE Is METERS 
| BEGINNING WITH SURFACE (PROGRAM 
SlepomernNo> MEreR SINCREMENTS 
DOWN TO 200 METERS) 
GRECULATE SOUND VELOCITY IN METERS 
PER SECOND FOR EACH 1 METER 
INCREMENT OF DEPTH AND STORE 
ANSWER IN TEMPORARY ARRAY 


RoiURN oe 
MENU 
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4, Inout Sea Surface Temperature And Average Gradient 


INPUT SALINITY IN 
PARTS PER THOUSAND FOR 
OPERATING AREA, ENTER A 
VALUE OF 35 IF SALINITY 
IS UNKNOWN 














INPUT LATITUDE IN WHOLE 


| 
DEGRZeo FOR OPERATING AREA 









TNEUTwoE. SURFACE 
TEMPERATURE IN 
VESGnooomenlolls 


INPUT AVERAGE TEMPERATURE 
GRADIENT IN DEGREES 
MSS IGS VSR bape | BON. 


me sLOWSS Ge OPP inh OF ENIERES) 
| (PROGRAM USES 200 METERS 
RomeOn fon Oot. ao 


CALCU DA ReaSOUND VELOCITY 
Pies Sonoran SECOND hOKR 
EACH 1 METER INCREMENT OP 
DEPEO AND STORE IN 
TEMPORARY ARRAY 


RETURN 
LO. MENU 


So 








Se 


Input Sound Ve locieyerrota le 












EINPG? SOUND VELOGCETY 
IN METERS] PER SEE€e@ND FOR 
INCREMENTS OF DEPTH IN METERS 
(PROGRAM USED 5 METER 
INCREMENTS) sslORE. VELGGITIES 
IN A TEMPORARY ARRAY 
(PROGRAM INTERPOLATES 1 METER 
INCREMENTS FROM 5 METER INPUTS) 
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6. Compute Range vs. Depth Solution And Store In Random 
Recess sie 


(5) 


INPUT HYDROPHONE/RECEIVER 
DEPTH IN EITHER 
FEET OR METERS 


| LOOP THROUGH EACH POSSIBLE 
|'SOURCE DEPTH, (THIS PROGRAM 
iUSES 10 METER INCREMENTS OF 
' DEPTH DOWN TO 200 METERS) 
COMPUTING PATH LENGTH 
DIFFERENCES FOR EACH 
INCREMENT OF HORIZONTAL 
RANGE, (THIS PROGRAM USES 
10 METER INCREMENTS OF 
HORIZONTAL RANGE OUT 
TO 1000 METERS) 
BY DOING THE FOLLOWING: 


I, SDEGERM ENE THE .UDP ER 
AND LOWER LIMITS FOR THE 
FOUR SEPARATE NUMERICAL 
INTEGRATIONS 


2. LOOP THROUGH EACH 
| POSSIBLE INITIAL RAY 
| ANGLE FROM EACH SOURCE : 
| DEPTE COMPUTING THE FOLLOWING: | 


RETRIEVE THE 
REQUIRED SOUND VELOCITIES 
FOR EACH INTEGRATION 


COMPUTE EACH 

UNIQUE SNELL CONSTANT 
GIVEN AN INITIAL RAY 
ANGLE AND SOUND VELOCITY 
DORgeeas SOURCE DEP TA 


NUMERICALLY SOLVE 

THE FOUR INTEGRALS DERIVED 
IN THE MATHEMATICAL 
SOLUTION DISCUSSION, 





cont. p. 
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————-~A CON. £rEom p. 6 


USING A ROMBERG IMPROVE- 
MENGE REE. 7 SON THE 
TRAPAZOIDAL NUMERICAL 
INTEGRATION TECHNIQUE 


IF INTEGRATION ERROR ON 
OCCURS WARN USER AND 
RETURN TO MENU ERROR 


STORE EACH INTEGRATION 
ANSWER IN THE APPROPRIATE 
TEMPORARY ARRAY 


CONTINUE TO LOOP THROUGH 
ALL POSSIBLE RAY ANGLES 
PRO aides SOURCE: (28s Obs 
PROGRAM USES ANGLES FROM 
3 DEGREES FROM THE VERT- 
ICAL TO 177 DEGREES FROM 
JHE WERTICAL, IN 3 DEGREE 
INCREMENTS ) 


USE LINEAR INTERPOLATION TO 
CALCU: THE £( DEP rERENCE) 


FOR 


EACH INCREMENT OF RANGE 


POR AYGivEN, DEPTE 


STORE ANSWERS IN RANDOM 
ACCESS FILE 


CONTINUE TO LOOP THROUGH 
DEPTH INCREMENTS UNTIL 


DONE 


RETURN 
TQ MENU 


(B 
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vy. Dasplayekinswer PisliemoOnsGr. Oz List.On PEinter 


DISPLAY OUTPUT MENU TO 


USER ALLOWING FOR SELECTION 
OMe 

TV. “GRay SGREEN 

Zs S2eINTER 





RETRIEVE HORIZONTAL RANGE VS. 
DEPTH ANSWER MATRIX FROM 
RANDOM ACCESS FILE AND 
OUTZ CG ® USERS, CHOICE 
FROM ABOVE 


t 









RETURN 
TO MENU 
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8. Exit Program And Return To Operating system 





WARN USER THAT THIS 
OPERATION WILL DESTROY 
ALL TEMPORARY FILES 
SUCH AS TEMPERATURE 
AND SOUND VELOCITY 
PROFILES 


ASK USER IF 
Enis Sat, NO RETURN 
DESIRED, YES OR NO TO MENU 


YES (B) 


EXIT TO OPERATING 
SEES OY 
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Wel 


ADD 


AFLAG 


ANS 


AOUT 


BFLAG 


BNDLOW 


BOUNDHI 
BOUNDLO 


CLAT 


DDEEP 


DEEP 


DEG. 


DELTA 


DHI 


DIRFRAC 


DIRRAY 


DEFINITION of COMPUTER PROGRAM VARIABLES 


Equation used to compute raypath length in 
meters for both the direct and reflected paths. 


Branching flag used in solving for answer matrix 
elements 


Answer obtained from numeric integration 


Answer matrix element output to CRT screen or 
printer 


Branching flag used in solving for answer matrix 
elements 


Lower limit for evaluating numeric integration 
bel 


Upper limit for evaluating numeric integration 
Direct path lower limit 


Latitude variation used in sound velocity 
equation 


First half of depth effect used in sound 
velocity equation 


Second half of depth effect used in sound 
velocity equation 


Average whole degree latitude of operating area 
used to calculate CLAT 


Line integral uniform divisions used in Romberg 
numeric integration routine 


Direct path horizontal range upper bound used in 
interpolating direct path length for desired 
range increment 


Direct path interpolation fraction 


Direct path interpolated path length 
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DLONG 


DLOW 


ert 


DSHORT 


FANSWER.FL1 


FINAL ( 


GRAD 
Heebner P 
peo P TH 
PoP D 
ua 


IFLAG 


el 
IN 
INCR 


INN 


JMP 


METER 


) 


Berce pai LenGgth upper bound used in 
interpolating path length for desired range 
increment 

Direct path hoirzontal range lower bound used in 
Pitereolatingudimeet path length for desired 
range increment 

Depth increment used for outputing answer matrix 
Direct path length lower bound used in 
interpolating path length for desired range 
increment 


Random access file used to store the horizontal 
range vs. depth answer matrix 


Answer array used to temporarily hold the 
interpolated values before storing them in 
random access file 

Average temperature gradilent input by user 
Hydrophone/receiver depth in feet 
Hydrophone/receiver depth in meters 

Sound velocity at hydrophone/receiver depth 


Loop increment 


Loop increment used to branch through one of 
four integrated equations 


Loop increment 

Loop increment 

Loop increment 

Loop increment 
Counter 

Array element pointer 
Counter 

Counter 


Array element pointer 
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MSTEP 
MSV 
MTRHI 
MTRLO 
MTRS 
MTRSTEP 
MULT 

N 
NANGL 
NMTR 
NN 
NSTEP 
NT 
NTRA 
NUM 
NUMB 
NUPLMT 


NX  ( 


PART 


PREGES 


RAD 


RAYDIR 


RAYREF 


Loop increment 

Array element pointer 
Array element pointer 
Array element pointer 
Array element pointer 
Array element pointer 
Random access file pointer 
Array element pointer 

Ray angle from source 
Horizontal range increment 
Counter 

Ray angle from source 
Array element pointer 
Array element pointer 

Main menu index selection 
Output menu index selection 
Incremented source depth 


Mmaeay of counters used during mumtegqration 
accuracy testing loop 


Bastitiloming factor used in numerical 
integration routine 


Incremented divisor used to partition depth 
integration limits 


Conversion of ray angle from degrees to radians 


Direct ray path length calculated from equation 
#12 page 25 


Reflected path length calculated from equation 
#13 page 25 
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REFFRAC 


REFLO 
REFRAY 
RHI 
RLONG 
RLOW 
RMTR 
RNGDIR 
RNGREF 


RSHORT 


SALIN 
SNDHI 
SNDLO 


SNDSPS 


SNELL 


SOUND 
Sst 
START 


SUM 


Reflected path interpolation fraction 


Reflected path lower limit. Negative value of 
boundlo 


Reflected path interpolated path length 
Reflected path horizontal range upper bound used 
in interpolating reflected path length for 
desired range 

Reflected path length upper bound used in 
interpolating path length for desired range 
increment 

Reflected path horizontal range lower bound used 
in interpolating reflected path length for 
desired range 


Depth increment for interpolation routine 


Direct path horizontal range array calculated 
from equation #7 page 21 


Reflected path horizontal range array calculated 
from equation #9 page 22 


Reflected path length lower bound used in 
interpolating path length 


Salinity value for operating area input by user 
Sound velocity at numeric integral upper bound 
Sound velocity at numeric integral lower bound 
Sound velocity profile array input directly by 
user or calculated from temperature, salinity, 


and latitude inputs 


Snell constant equation result calculated from 
initial ray angle and depth of source 


Sound velocity at incremented source depth 
Sea surface temperature 
Initial estimate of integral answer 


Equation used to compute horizontal range in 
meters for both the direct and reflected paths 
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SVP 
oor 


BEST ( 


TGTDPTH 


TOL 


mer 


TTEMP 


WALK 


) 


Input sound velocity profile by user 
Degrees celsius temperature divided by 10 


Successive numeric integration answer array used 
to check accuracy of answer 


Incremented source depth 


Accuracy desired for numerical integration 
routine answer 


Summation of partitioned numerical ilntegration 
values 


Input temperature used in calculating sound 
velocity profile 


Integration routine smoothing factor 


Incremented depth value for each partition 
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APPENDIX A 


Urick [Ref. 1: p. 124] briefly discusses the intensity of sound 
reaching a receiver via a direct path and a surface reflected path. The 
following derives the ratio of the total intensity of the signal received 
at a point when multipath Lloyd's mirror interference is present, to the 


total intensity resulting only from the direct path. 


RANGE 
METERS 





R 


Zz DIRECT ioe a> 
S 
DEPTH Intensity Ratio = 5 » 
METERS + 0 0 
MUST SHOW: 
Ry \2 R 
<I> 2iieyD D 
sa = ltypltsj] - af?) cos(wt ) 
as (re) Rp 


WHERE: <I>, the time averaged intensity and I, the instantaneous 


intensity are defined as follows: 


. 
Q) <b arn) f I dt T = 1 period 
wa 
0 aig 
(2) (Po(r) + Pa(r))° ~ 
I = sea Pa(r) and Pa(r) are both functions of time 
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GIVEN: 


= P P 


(1) Peecervep * Porrect * PrerLectep 


Po uP 9 
R— cos wt - g— cos [a(t-7)] 
D R 


where: wu = surface pressure reflection coefficient 


of 
bs 


time delay caused by the difference in the two path 


lengths 
p2 PG 
(2) Thy Ty emi 
2ocRy 


where: P = RMS pressure amplitude of the direct path ray 
NOTE: It is assumed that the transmission loss is due only to spherical 


Spreading and that absorption in the medium is negligible. 





SOLUTION: 
1) Oe UG 
POIR = . ( wt) Peer = R, cos w(t-r)] 
PRECEIVED ~ "prrecT * PREFLECTED 


PRECETVED = roc cos(ut Re cos (at-2) 


2 
2 ?P 2 2 
ai 0 | eons u 2 _ 
I = Je 72 + m. cos (uw(t-t)) 
D R 


- eu wt [w (t-t) ] 
RoR cos (wt) cos [w(t-rt) 





Sal 





SINCE: png [1 dt 


AND w= anf = & poe 
70) 


On/w 


2 
cos¢(w(t-t)) dt 


5 
RR 





a> = ia cosé(wt) dt + 
pC] p 
0 


2n/w 
-_—_ cos(wt) cos(w(t-1t)) dt 


Solving the three integrals separately 


2n/w 2n/w 


-- cos@(ut) dt = +, § + i sin(2 4 
RD : 0 


S2 





ale 2n/w 








2 2 = \fee fl 
3 cos (w(t-t)) dt = 255 pia es Sin(2 w(t-1)) 
R R A 
R Jp R 0 
nr, l 
= 5 a + i (sin(4n - 2y7)-sin(-2y7)) 
RR 
=| uam 
wre 
2n/w 
a cos(wt) cos(w(t-r)) dt 
RD 
0 
2m/ w 
= a cos“ (wt) cos(wt)+sin(wt) cos(wt) sin(wr) dt 
RD 
0 
on /w 2n/y 
= Sy, caster) cos“ (wt) dt + oy Sinkut) Sin(wt) cos(yt) dt 
RD RD 
0 0 
on /w ony 
= Su costut) (£4). sin(auty | + 2 Simlur) (sin?(,.t)) 
Rp Rp *® 
0 ) 
= 2u_cos (wt) (ny , 2u sin(wt) (9) 
RoR Ww ReRp 


2 On 
ORR cos(ut ) 


a8 





Assembling the answers to the three integrals yields: 




















“I> = oP sun 2m : 
Pape | 2 2 ~ wROR, COSCwt) 
wR wR RD 
D R 
2 
p 2 
3 US a aaa! 
a =. 2 * 2 ROR COS(wr ) 
yO 
pé Re 2uRy 
<I> = 1 + ~ COS( wr) 
20cR° po 
aa) R 
NOTE: 
Po 
<I> = 7 
20cRp 
ak cur, 
<I> = <I,> | 1 +—— - cos(wrt) 
0 we Ro 
R 


Finally 
oe =1l+ ue 2) - aul cos(wr) 
0 R R 
This result shows the dependence of the intensity at the receiver on 
the difference in path lengths RO IRECT and ReeriEcTED and the time delay, 
r, caused by this difference. This intensity ratio is a minimum (corre- 
sponding to a null on the Lloyd's mirror interference pattern) when 


cos wrt = 1 or wrt = integer multiples of 2m. 
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COMPUTER OULPU a 


A sample of the actual computer output is listed below. 


Data values are frequency measurements Af (Hertz) between 


Successive nulls on the Lloyd's mirror pattern (see Figure 


TI-1). 


Given a particular source depth and a Af, the source range 


can be determined from the corresponding output column. 


The 


source depth can be determined from the appropriate output row 


if the horizontal range to the source from the receiver and 


the Af are known. 


Group of Ranges From 10 to 70 Meters 


Depth of 
S056 


Depth of 
iS 


Depth of 
aS 


Depth of 
2052 


Depth of 
2B D 


10 Meters 


Oe 8 i dee 


20 Meters 
BO. 2 


30 Meters 
lien 5 


40 Meters 
SAS 


50 Meters 
30.6 


BOR 


43.8 


36.8 


205 


3167550 


alias 


SU 0 


41.7 


5629 


313) 


LESS Bes: 


BS 2 ent 


sien © 


46.9 


40.8 


184.4 


94.8 


66.1 


Sas) 


44.9 


Pa Oe 


Oye 3 


14.3 


Seg & 
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COMPUTER PROGRAM 


1O) DERINT 1,J0,K,L,M,N 

20 DIM NX(16),TEST(136),SNDSPD( 201),RNGDIR(180) 

30 DIM RNGREF(180),RAYDIR(180),RAYREF(180),F INAL(1 01) 

40 TOm= .05 

50 OPEN "R",#1,"FANSWER.FL1", 4 

60 FIELD #41, 4 AS F1$ 

100 PRINT CHR$(26): PRINT "SELECT OPERATION": PRINT 

110 PRINT 1,"INITIALIZE ANSWER FILE,DO AT THE BEGINNING OF 
EACH NEW BT INPUT" 

120 PRINT 2,"INPUT BT INFO IN 5 METER INCREMENTS" 

130 PRINT 3,"INPUT SEA SURFACE TEMP AND GRADIENT FOR FIRST 
200 METERS" 

moore PROENT 4,°ENPUT SOUND VELOCITY PROFILE IN 5 METER 
INCREMENTS" 

140 PRINT 5,"CALCULATE LLOYD'S MIRROR HORIZONTAL RANGE .VS. 
DEPTH CHART" 

150 PRINT 6,"OUTPUT CHART TO PRINTER" 

160 PRINT 7,"EXIT PROGRAM" 

170 PRINT:PRINT:INPUT "ENTER THE NUMBER OF THE OPERATION 
DESIRED " NUM 

180 IF (NUM<1) OR (NUM>7) THEN PRINT “IMPROPER INPUT":GOTO 
100 

190 ON NUM GOSUB 1000,2000,5000,5500, 6000, 9600,10000 

200 GOTO 100 

1000 REM INITIALIZE ANSWER FILE 

1010 PRINT CHR$(26):PRINT "ARE YOU SURE YOU WANT TO ZERO OUT 
THE ANSWER FILE" 

1030 INPUT "ENTER Y FOR YES, N FOR NO YES 

1040 IF A$<>"Y" THEN RETURN 

1050 CLOSE #1 

1060 KILL "FANSWER.FL1" 

1070 OPEN "R",#1,"FANSWER.FL1", 4 

1080 FIELD #1, 4 AS F1$ 

1090 RETURN 

2000 REM INPUT BT INFO IN 5 METER INCREMENTS, THEN CALCULATE 
SOUND VELOCITY 

2010 REM USING ALAN R. COPPENS EQUATION PRESENTED IN J. 
ACOUS. SOC. OF AMER. 

2015 REM VOLUME 69, PAGE 862, 1981 AND PUBLISHED IN 
FUNDAMENTALS OF ACOUSTICS 

2016 REM THIRD EDITION, BY KINSLER FREY COPPENS AND 
SANDERS 

2020 GOSUB 3010: GOSUB 4000:PRINT CHR3$( 26) 

2030 PRINT “ENTER BT INFO IN 5 METER INCREMENTS STARTING 
WITH THE SURFACE" 

2040 FOR I = 0 TO 200 STEP 5 

2050 PRINT:PRINT I “METERS":INPUT "TEMP IN CELSIUS = 
" TTEMP:PRINT 

2060 PRINT "YOU INPUT “TTEMP" DEGREES CELSIUS": 
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: INPUT "IS THIS CORRECT ? ENTER Y FOR YES N FOR NO 
7 BS 
2070 IF B$<>"Y" THEN PRINT "REENTER TEMP":GOTO 2050 
2075 TEMP = TTEMP / 10! 
2080 METER = I 
2100 DEEP = (16.23 +.253*TEMP)* CLAT * CSNG(I)/1000! +(.213- 
.1*TEMP) * 
(CLAT * (CSNG(I)/1000!))?2 
2105 DDEEP = (.016 + .0002 * (SALIN-35)) * (SALIN- 
35) *TEMP*CLAT*CSNG(1I)/1000! 
2110 SNDSPD(METER)=1449.05 + 45.7 * TEMP - 5.21 * (TEMP?2) + 
323 * (TEMP73) + 
(1.333-.126*TEMP+.009*(TEMP?2))*(SALIN- 
35)+(DEEP+DDEEP) 
Zin «6©Eet = 0 GOTO 2120 
2112 SLOW = SNDSPD( I - 5 ) 
2113 SHI = SNDSPD( I ) 
maBOR J =) 1 TO 4 
2115 SNDSPD( I - 5 +J ) = SLOW +(J9/5) * ( SHI - SLOW ) 
216 NEXT J 
2120 NEXT I 
2130 RETURN 
3000 REM ENTRY OF OPAREA LATITUDE 
3010 PRINT CHR$(26):PRINT "ENTER AVERAGE DEGREES LATITUDE OF 
OPAREA" 
3020 INPUT "IN WHOLE DEGREES " DEG 
3030 PRINT:PRINT "YOU ENTERED "DEG" DEGREES, IS THIS CORRECT 
att 
3040 INPUT "ENTER Y FOR YES N FOR NO HSEGS 
3050 IF cc$<>"yY" THEN PRINT "REENTER LATITUDE ":GOTO 3010 
3060 CLAT = 1 - .0026 * (COS( DEG * 3.14159 / 90!)) 
3070 RETURN 
4000 REM ENTRY OF SALINITY VALUE FOR SOUND SPEED 
COMPUTATIONS 
4010 PRINT CHR$(26) 
4020 PRINT "ENTER THE SALINITY VALUE AS PARTS PER THOUSAND 
FOR THE OPAREA" 
4030 INPUT "IF UNKNOWN, ENTER A VALUE OF 35, SALINITY = 
" SALIN 
4040 PRINT:PRINT "THE SALINITY LEVEL IS "SALIN" PARTS PER 
THOUSAND" 
4050 INPUT "IS THIS CORRECT ? ENTER Y FOR YES N FOR NO 
" BBS 
4060 IF BB$ <> "Y" GOTO 4010 
4070 RETURN 
5000 REM ENTER SEA SURFACE TEMP AND AVERAGE GRADIENT FOR 
FIRST 200 METERS 
5010 PRINT CHR§$( 26) 
5020 INPUT "ENTER THE SEA SURFACE TEMPERATURE IN DEGREES 
CELSIUS " SST 
5030 PRINT:PRINT "YOU ENTERED "SST" DEGREES CELSIUS FOR THE 
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SEA SURFACE TEMP" 
5040 INPUT "IS THIS CORRECT ? ENTER Y FOR YES N FOR NO 
ae © 
5050 IF c$ <> "Y" THEN PRINT "REENTER SEA SURFACE TEMP":GOTO 
5020 
5060 PRINT 
5070 PRINT “ENTER THE TEMPERATURE GRADIENT (DEG 
CELSIUS/METER)" 
5075 INPUT "FOR THE FIRST 200 METERS " GRAD: PRINT 
5080 PRINT "YOU ENTERED A TEMPERATURE GRADIENT OF "GRAD" 
DEGREES CELSIUS/METER" 
5090 INPUT “IS THIS CORRECT ? ENTER Y FOR YES N FOR NO 
,D$ 
5100 IF D$ <> 'Y¥Y" THEN PRINT "REENTER GRADIENT":GOTO 5070 
5110 GOSUB 3010 
5120 GOSUB 4000 
5125 PRINT CHR$(26):PRINT "PLEASE WAIT, CALCULATING SOUND 
VELOCITY VALUES" 
5130 FOR I = 0 TO 200 
5140 TEMP = (SST + (CSNG(I) * GRAD)) / 10! 
5150 DEEP = (16.23 + .253*TEMP)* CLAT * CSNG(I) /1000! + 
(.213-.1*TEMP) 
wo LAr * (CSNG(I) /1000!: ))72 
5155 DDEEP = (.016 + .0002 * (SALIN - 35)) * (SALIN-35) * 
TEMP * CLAT * 

CSNG(I) / 1000! 


5170 SNDSPD(I) = 1449.05 + 45.7 * TEMP - 5.21 * (TEMP?2) + 
os (PRP? 3) 
+#(1.333 - .126 * TEMP + ,.009*TEMP?2)*(SALIN - 


35)+DEEP+DDEEP 

5180 NEXT I 

5190 RETURN 

5510 REM ENTRY OF SOUND VELOCITY PROFILE 

5525 PRINT CHR$(26) 

5530 PRINT “ENTER SOUND VELOCITY PROFILE IN 5 METER 
INCREMENTS " 

5540 PRINT "STARTING WITH THE SURFACE SOUND VELOCITY IN 


METERS/SECOND" 

5550 FOR II = 0 TO 200 STEP 5 

§560 PRINT: PRINT It “METERS“" :INPUT " SOUND 
VELOCITY(METERS/SEC)= ",SVP 


5570 PRINT: PRINT "YOU INPUT "SVP" METERS/SEC" 
5580 INPUT "IS THIS CORRECT ? ENTER Y FOR YES N FOR NO 
"eSS$ 
rf 
5590 IF sSS$<>"Y" THEN PRINT "REENTER SOUND VELOCITY": GOTO 
5560 
5610 SNDSPD(II) = SVP 
5611 IF II = 0 GOTO 5620 
5612 SLOW = SNDSPD( II - 5 ) 
5613 SHI = SNDSPD( II ) 
5614 FOR JJ = 1 TO 4 
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5615 SNDSPD( II - 5 + JJ ) = SLOW + ( JJ / 5 ) * ( SHI - 
SLOW ) 

5616 NEXT JJ 

5620 NEXT II 

5630 RETURN 

6000 REM CALCULATE LLOYD'S MIRROR HORIZONTAL RANGE .VS. 
DEPTH CHART 

6005 PRINT:PRINT:PRINT "WILL HYDROPHONE DEPTH BE INPUT IN 
FEET OR METERS ?" 

6006 INPUT “ENTER F FOR FEET OR M FOR METERS " AHYD$ 
6007 IF AHYD$ = "F" GOTO 6010 ELSE IF AHYD$ = "M" GOTO 6052 
6008 PRINT: PRINT: PRINT "INPUT ERROR" :GOTO 6005 

6010 PRINT CHR$(26):INPUT "ENTER HYDROPHONE DEPTH IN FEET 
" HYDEEP 

6020 PRINT:PRINT "YOU ENTERED A HYDROPHONE DEPTH OF "HYDEEP" 
FEET" 

6030 INPUT "IS THIS CORRECT ? ENTER Y FOR YES N FOR NO 
" AAS 

6040 IF AA$ <> "Y¥" GOTO 6010 

6050 HYDPTH = (HYDEEP * .3048) 

6051 GOTO 6060 

6052 PRINT CHRS(26):INPUT "ENTER HYDROPHONE DEPTH IN METERS 
" HYDPTH 

6053 PRINT:PRINT "YOU ENTERED A HYDROPHONE DEPTH OF 
"HYDPTH" METERS" 

6054 INPUT " IS THIS CORRECT ? ENTER Y FOR YES N FOR NO 
" BHY$ 

6055 IF BHY$ <> "Y" GOTO 6052 

6060 HYSPD = SNDSPD( HYDPTH ) 

6170 FOR NUPLMT = 10 TO 200 STEP 10 

6175 PRINT "NUPLMT=" NUPLMT 

6180 TGTDPTH = CSNG( NUPLMT ) 

6185 IF HYDPTH = TGTDPTH THEN HYDPTH = HYDPTH - 3 

6190 IF TGTDPTH > HYDPTH GOTO 6230 

6200 BOUNDLO = TGTDPTH 

6210 BOUNDHI = HYDPTH 

6220 GOTO 6250 

6230 BOUNDLO = HYDPTH 

6240 BOUNDHI = TGTDPTH 

6250 REFLO = -(BOUNDLO) 

6360 FOR NSTEP = 3 TO 177 STEP 3 

6365 PRINT "NSTEP=" NSTEP 

6370 RAD = CSNG(NSTEP) * 3.14159 / 180! 

6380 MTRSTEP = TGTDPTH 


6390 SOUND = SNDSPD(MTRSTEP) 
6400 MTRHI = BOUNDHI 

6410 SNDHI = SNDSPD(MTRHI) 

6420 MTRLO = BOUNDLO 

6430 SNDLO = SNDSPD(MTRLO) 

6440 SNELL = SIN( RAD ) / SOUND 


6480 FOR IFLAG = 1 TO 4 
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7000 REM ROMBERG NUMERICAL INTEGRATION SUBROUTINE 
70g SUART = 0! 

7020 PIECES = 1! 

Oe NX(1) = 1 

7040 IF (IFLAG=2) OR (IFLAG=4) THEN BNDLOW=REFLO ELSE 
BNDLOW=BOUNDLO 

7041 IF (IFLAG = 2) OR (IFLAG = 4) GOTO 7045 

7042 IF (TGTDPTH > HYDPTH) AND (NSTEP > 90) GOTO 8105 
7043 IF (TGTDPTH < HYDPTH) AND (NSTEP < 90) GOTO 8105 
7044 GOTO 7050 

7045 IF NSTEP > 90 GOTO 8105 

7050 DELTA = (BOUNDHI - BNDLOW) / PIECES 

7051 DENOM1 = (SNELL)?2 * (SNDLO) ?2 

7052 DENOM2 = (SNELL)?2 * (SNDHI) 72 


Poe ele sPENOMI = 1! THEN DENOM1 = .999 :1F DENOM2 = 1! THEN 
DENOM2 = .999 
7054 IF ( DENOM1 > 1! ) OR ( DENOM2 > 1! ) GOTO 8105 


7060 IF (IFLAG=1) OR (IFLAG=2) GOTO 7085 

7065 FOR IN = 1 TO 2 

7068 IF IN = 1 THEN DENOM = DENOM1 ELSE DENOM = DENOM2 
7070 START = START + ((1!)/ SQR(1! - DENOM)) 

7075 NEXT IN 

7080 GOTO 7100 

7085 FOR INN = 1 TO 2 

7088 IF INN = 1 THEN DENOM = DENOM1 ELSE DENOM = DENOM2 
7089 IF INN = 1 THEN SPEED = SNDLO ELSE SPEED = SNDHI 
7090 START = START + ((SNELL * SPEED) / SQR(1! —- DENOM)) 
7095 NEXT INN 

7100 TOT = START / 2! 

Tare GPEST( 1!) = DELTA * START / 2! 

7120 N = 1 

7130 NN = 2 

7140 Ne= N + 1 

7150 WALK = 4! 

7160 NX(N) = NN 

7170 MSTEP = CINT( PIECES ) 

7180 PIECES = PIECES * 2! 

7190 DELTA = (BOUNDHI - BNDLOW) / PIECES 


7500 REM COMPUTE TRAPEZOIDAL SUM 
Zoe BOR ENCR = 1 TO MSTEP 
Voc) PART = CSNG(INCR) * 2! - 1! 


7530 X = BNDLOW + (DELTA * PART) 

7540 MTRS = CINT ( ABS (X) ) 

7545 IF MTRS > 200 THEN MTRS = 200 

7550 SPEED = SNDSPD( MTRS ) 

7552 DENOM3 = (SNELL)?2 * (SPEED) ?2 

7554 IF DENOM3 = 1! THEN DENOM3 = .999 

7555 IF DENOM3 > 1! GOTO 8105 

7570 IF (IFLAG=1) OR (IFLAG=2) GOTO 7600 
7580 TOT = TOT + ( (1!) / SQR(1! - DENOM3)) 
7590 GOTO 7610 
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7600 
7610 
7620 
7630 
76 40 
7650 
7660 
7670 
7680 
6 FO 
7700 
7705 
8105 
i 10 
7720 
7730 


TOT = TOT + ( (SNELL * SPEED) / SQR(1! - DENOM3)) 
NEXT INCR 
TEST( NN ) = 
="N - 7 
NTRA = NX( 
FOR M = 1 
J = NN +M 
NT = NTRA + M - 1 
TEST(J) = ( TEST(J-1) 
WALK = WALK * 4! 
NEXT M 

IF (IFLAG = 


TOT * DELTA 


K ) 
TO K 


* WALK - TEST(NT)) / (WALK-1!) 


1 OR IFLAG = 2) AND (TEST(J) >» 3000!) GOTO 


IF N < 5 GOTO 7760 
mre ort. NN + 1) = 0! GOTO 7750 
IF ABS(TEST(NTRA+1)-TEST(NN + 1)) 


1)*TOL)) GOTO 8000 


7740 
GOTO 
7750 


8110 
6115 


IF ABS( TEST(NN-1) - 
8000 

Lr N > 15 GOTO 9500 
Nie= J + 1 

GOTO 7140 

ANS = TEST(J) 

ON IFLAG GOTO 8040, 8060, 8080, 8100 
RNGDIR(NSTEP) ANS 

GOTO 8105 
RNGREF (NSTEP ) 
GOTO 8105 
RAYDIR(NSTEP) 
GOTO 8105 
RAYREF (NSTEP ) 
NEXT IFLAG 
NEXT NSTEP 
REM USE LINEAR INTERPOLATION TO CALCULATE ANSWER CHART 


TEST(J)) <= (ABS( TEST(J) * TOL)) 


ANS 


ANS 


ANS 


VALUES IN RANGE 


ory 
8120 
8130 
€135 
8140 
on 50 
8160 
8161 
8162 
o163 
8164 
8165 
8166 
8167 
8168 
8169 


REM AND DEPTH INCREMENTS OF 10 METERS 
FOR NMTR = 1 TO 100 

RMTR = CSNG( NMTR ) * 10! 

DIRRAY = 0! : REFRAY = 0! 

AFLAG = Q! 

BFLAG = 0! 


IF TGTDPTH < HYDPTH GOTO 8179 


FOR NANGL = 3 TO 87 STEP 3 

IF RNGDIR( NANGL ) <= RMTR GOTO 8171 

DLOW = RNGDIR( NANGL - 3 ):DHI = RNGDIR( NANGL ) 
IF DLOW = DHI THEN DIRFRAC = 0 :GOTO 8166 
DIRFRAC = (RMTR - DLOW) / (DHI - DLOW) 


IF DLOW > 0 GOTO 8168 


DSHORT = BOUNDHI - BOUNDLO :GOTO 8169 
DSHORT = RAYDIR( NANGL - 3 ) 
DLONG = RAYDIR( NANGL ) 
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Bil 70 
8210 
ot 
ELSE 
8172 
e173 
8174 
e175 


DERRAY = DSHORT + ( DLONG — DSHORT ) * DIRFRAC : GOTO 


IF ( NANGL = 87 ) AND ( RNGDIR( 87 ) > 0 ) GOTO 8172 
GOTO 8177 

DLOW = RNGDIR( 87 ) 

IF DLOW = RMTR THEN DIRRAY = RAYDIR( 87 ) : GOTO 8210 
DIRFRAC ( RMTR - DLOW ) / ( 1000 - DLOW ) 

DSHORT = RAYDIR( NANGL ) : CLONG = SQR( 100072 + 


(BOUNDHI - BOUNDLO)?72 ) 


8176 
Sz 100 
O77 
8178 
St 9) 
8180 
8181 
ELSE 
8182 
S11 OS 
8184 


DIRRAY = DSHORT + ( CLONG - DSHORT ) * DIRFRAC : GOTO 


NEXT NANGL 

GOTO 8410 

FOR MANGL = 93 TO 177 STEP 3 

IF RMTR >= RNGDIR( MANGL ) GOTO 8181 ELSE GOTO 8195 

IF ( MANGL = 93 ) AND ( RNGDIR( 93 ) > O ) GOTO 8182 
GOTO 8186 

DLOW = RNGDIR( 93 ) 

IF DLOW = RMTR THEN DIRRAY = RAYDIR( 93 ) : GOTO 8210 
DIRFRAC = ( RMTR - DLOW ) / (1000 - DLOW ) : DSHORT = 


RAYDIR( 93 ): 


8185 
8210 
8186 
8187 
GOTO 
8188 
8 139 
8190 
8191 
ag 
8194 
8210 
o.95 
8196 
81 Oe 
8198 
8210 
Bagg 
8200 
8205 
8210 
e211 
SIZ 2 
e713 
8214 
Be Wes, 
8216 
Syyayll 7 
e2rg 


BLONG = SOR ( 100072 + ( BOUNDHI - BOUNDLO )?72 ) 
DIRRAY= DSHORT + ( BLONG - DSHORT ) * DIRFRAC : GOTO 


DLOW = RNGDIR( MANGL ) :DHI = RNGDIR( MANGL - 3 ) 
IF ( DLOW <= RMTR ) AND ( DHI >= RMTR ) GOTO 8188 ELSE 


IF DLOW = DHI THEN DIRFRAC = 0 : GOTO 8190 

DIRFRAC = ( RMTR - DLOW ) / ( DHI - DLOW ) 

IF DLOW > 0 GOTO 8192 

DSHORT = BOUNDHI - BOUNDLO : GOTO 8193 

DSHORT = RAYDIR( MANGL ): DLONG = RAYDIR( MANGL - 3) 
DIRRAY = DSHORT + ( DLONG - DSHORT ) * DIRFRAC : GOTO 


IF ( MANGL <> 177 ) GOTO 8199 

DLOW = 0 : DHI = RNGDIR( 177 ) : DIRFRAC = RMTR / DHI 
DSHORT = BOUNDHI - BOUNDLO : DLONG = RAYDIR( 177 ) 
DIRRAY = DSHORT + ( DLONG - DSHORT ) * DIRFRAC : GOTO 


NEXT MANGL 

GOTO 8410 

IF DHI = DLOW THEN DIRFRAC = 0 : GOTO 8211 

FOR NREF = 3 TO 87 STEP 3 

IF RNGREF( NREF ) <= RMTR GOTO 8220 

RLOW = RNGREF( NREF - 3 ) : RHI = RNGREF( NREF ) 
IF RLOW = RHI THEN REFFRAC = 0 : GOTO 8215 
REFFRAC = ( RMTR - RLOW ) / ( RHI - RLOW ) 

iE Rowe 0 GOTO 8217 

RSHORT = BOUNDHI + BOUNDLO : GOTO 8218 

RSHORT = RAYREF( NREF - 3 ): RLONG = RAYREF(NREF) 
REFRAY = RSHORT + ( RLONG - RSHORT ) * REFFRAC :GOTO 
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8400 

8220 IF ( NREF = 87 ) AND (RNGREF(87) > O ) GOTO 8221 ELSE 
GOTO 8226 

8221 RLOW = RNGREF( 87 ) 

8222 IF RLOW RMTR THEN REFRAY = RAYREF( 87 ): GOTO 8400 
8223 REFFRAC ( RMTR - RLOW ) / ( 1000 - RLOW ) 

8224 RSHORT = RAYREF( NREF ) : ALONG = SOR ( 100072 + ({ 
BOUNDHI + BOUNDLO )?2 ) 

8225 REFRAY = RSHORT + ( ALONG - RSHORT ) * REFFRAC : GOTO 
8400 

8226 NEXT NREF 

8227 GOTO 8410 

8370 REM OUTPUT FINAL ANSWERS FOR EACH DEPTH AND RANGE TO 
DISK FILE 

8385 LPRINT "HYSPD=" HYSPD 


, i 


8400 FINAL( NMTR ) = HYSPD / ( REFRAY - DIRRAY ) 
8410 NEXT NMTR 
8420 MULT = (TGTDPTH + 10) - 1 


8430 FOR JMP = 17T0 100 
8440 LSET F1$ = MKS$(FINAL(JMP)) 
8450 RECNUMS = JMP + (MULT ¥* 100) 
8460 PUT #1, RECNUM$% 
8470 NEXT JMP 
8480 FOR NUL = O TO 100 
8490 FINAL(NUL) = 0! 
8500 NEXT NUL 
8510 FOR NUL = 0 TO 135 
8520 TEST(NUL) 0! 
8530 NEXT NUL 
8540 FOR NUL = O TO 15 
So SOmNAxt( 15) = 0! 
8560 NEXT NUL 
8570 FOR NUL = O TO 179 
85 80RNGDIR(NUL)=0 !:RNGREF(NUL)=0:RAYDIR(NUL)=0: 
RAYREF (NUL) =0 
8620 NEXT NUL 
9000 NEXT NUPLMT 
9010 RETURN 
9500 PRINT: PRINT:PRINT "INPUT DATA HAS CAUSED AN ERROR" 
9510 GOTO 100 
9600 REM PRINT OUT GRAPH OF HORIZONTAL RANGE .vs. DEPTH 
9610 PRINT CHR3$(26) 
9620 PRINT:PRINT:PRINT " SELECT OUTPUT DEVICE ": PRINT 
9630 PRINT 1," VIDIO SCREEN ": PRINT 
9640 PRINT 2," PRINTER ": PRINT 
9650 INPUT " ENTER THE NUMBER FOR OUTPUT DEVICE DESIRED 
" NUMB 
9660 IF ( NUMB < 1 OR NUMB > 2 ) THEN PRINT "IMPROPER INPUT 
"sGOTO 9620 
9665 PRINT:PRINT 
9670 PRINT " THE OUTPUT MATRIX IS PRESENTED IN 10 METER 


it) 
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INCREMENTS FOR BOTH " 
9680 PRINT " HORIZONTAL RANGES AND FOR DEPTHS OF 10 METERS 
TO 200 METERS, " 
9690 PRINT " HORIZONTAL RANGES FROM THE RECEIVER TO THE 
FARGET , IN ” 
9700 PRINT " GROUPS OF 100 METERS, INCREASE FROM LEFT TO 
RIGHT, 
9710 PRINT “ ie. 10 TO 100 METER RANGES WILL BE DISPLAYED 
PRS. KOLLOWED 
9720 PRINT " BY 110 TO 200 METERS AND SO ON UP TO 1000 
METERS." 
9730 PRINT " DEPTHS OF 0 TO 200 METERS WILL BE DISPLAYED 
FOR EACH GROUP OF " 
9740 PRINT " RANGES, WITH THE FIRST ROW REPRESENTING 10 
METERS, THE " 
9750 PRINT " NEXT ROW 20 METERS, ON DOWN THE PAGE UNTIL ROW 
20 IS 
9760 PRINT “" PRINTED FOR A DEPTH OF 200 METERS.":PRINT: PRINT 
9770 FOR I = 0 TO 9 
9780 THIS = 10 + 100 * I 
9790 THAT = 100 + 100 * [I 
9800 IF NUMB = 1 THEN PRINT: 
PRINT "GROUP OF RANGES FROM “THIS" TO 
"THAT"METERS":GOTO 9820 
9805 LPRINT:LPRINT: 
9810 LPRINT:LPRINT "GROUP OF RANGES FROM "THIS" TO "THAT" 
METERS" 
9820 FOR J = 0 TO 19 
9830 DPH = 10 + 10 * J 
9840 IF NUMB = 1 THEN PRINT:PRINT "DEPTH OF "DPH" 
METERS":GOTO 9860 
9850 LPRINT:LPRINT "DEPTH OF "DPH" METERS" 
9860 FOR K = 1 TO 10 
9870 RECNUMS = K+100* J +10* I 
9880 GET #41,RECNUM$% 
9890 AOUT = CVS( F1$ ) 
9900 IF NUMB = 2 GOTO 9930 
9910 PRINT USING "####.# '"sAOUT; 
9920 GOTO 9940 
9930 LPRINT USING "####.# '":AOUT; 
9940 NEXT K 
9950 NEXT J 
9960 NEXT I 
9965 LPRINT 
9970 RETURN 
10000 REM CLOSE ALL FILES AND EXIT PROGRAM 
10010 PRINT CHR$(26):PRINT "DO YOU REALLY WANT TO QUIT" 
10020 INPUT "ENTER A Y FOR YES, AN N FOR NO Zo 
10030 IF Z2$ = "Y'" THEN END ELSE RETURN 
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